\documentclass{biophys}
\usepackage{dcolumn,color,helvet,times}
\usepackage{amsmath,lastpage,bm,textcomp}%,mathtime}
\usepackage{multicol}
\usepackage{amssymb,amsthm}
%\usepackage{epsfig}
%\usepackage{showkeys}
\usepackage{diagrams}
\usepackage{amscd}
\usepackage{chemarrow}
\usepackage{amsmath,amsfonts} % for example for align command
\usepackage{booktabs} % for much better looking tables
\usepackage{array} % for better arrays (eg matrices) in maths
\usepackage{paralist} % very flexible & customisable lists (eg. enumerate/itemize, etc.)
\usepackage{verbatim} % adds environment for commenting out blocks of text & for better verbatim
%\usepackage{subfig} % make it possible to include more than one captioned figure/table in a single float
% These packages are all incorporated in the memoir class to one degree or another...

\usepackage[numbers]{natbib}
\bibliographystyle{chicagoa}
\usepackage{graphicx} % support the \includegraphics command and options
%\usepackage[parfill]{parskip} % Activate to begin paragraphs with an empty line rather than an indent

\jno{} %journal number
\gridframe{N}%option for grid around the text "Y" or "N"
\cropmark{N}%option for cropmark around the text "Y" or "N"
\doi{}% DOI number in the copyright line

\makeatletter
\newenvironment{tablehere}
  {\def\@captype{table}}
  {}

\newenvironment{figurehere}
  {\def\@captype{figure}}
  {}
\makeatother


\newcommand{\md}{d\kern-0.035cm\char39\kern-0.03cm}
\newcommand{\mL}{L\kern-0.035cm\char39\kern-0.03cm}
\newcommand{\eq}[1]{(\ref{#1})}
\newtheorem{Theorem}{Theorem}
\newtheorem{Definition}{Definition}


\newcommand{\la}{\lambda}
\newcommand{\eps}{\varepsilon}
\newcommand{\el}{\ell}
\newcommand{\ksi}{\xi}

\newcommand{\sarr}{\rightarrow}
\newcommand{\darr}{\rightleftarrows}
\newcommand{\dsum}{\displaystyle \sum}

\newcommand{\CC}{\mathbb{C}}
\newcommand{\C}{\CC}
\newcommand{\RR}{\mathbb{R}}

\newcommand{\dG}{\triangle G}

\setcounter{page}{1} %first page number

\title{Mathematical model of alternative mechanism of telomere length maintenance in yeast mtDNA}
\author{Richard Koll\'ar*, Katar{\'\i}na Bo{\md}ov\'a*, Jozef Nosek{\authdagger}, 
and {\mL}ubom{\'\i}r Tom\'a{\v s}ka{\authdagger}{\authdagger}}

\address{*Department of Applied Mathematics and Statistics, Faculty of Mathematics, Physics, 
and Informatics, Mlynsk{\'a} dolina, 842 48 Bratislava, Slovakia, {\addrdagger}
Department of Biophysics and {\addrdagger}{\addrdagger}Department of Genetics, 
Faculty of Natural Sciences, Mlynsk{\'a} dolina, 842 15 Bratislava, Slovakia} 
\date{\today} % Activate to display a given date or no date (if empty),
         % otherwise the current date is printed  
\markboth{Biophysical Journal: Biophysical Letters}{Biophysical Journal: Biophysical Letters} %for running head


\begin{document}
%\onecolumn
\maketitle

\pagestyle{headings}

%Abstract environment needs 3 arguments. They are
%1. The abstract
%2. Received date
%3. Address, email

\begin{abstract}
{Telomeres, non-coding terminal structures of DNA strands, consist of repetitive long tandem repeats of a specific length. 
Absence of the enzyme telomerase in certain cellular structures requires an alternative telomerase-independent 
pathway for telomeric sequence length regulation. Besides linear telomeres other configurations such as telomeric 
circles and telomeric loops have been experimentally observed. 
We propose a mathematical model that captures biophysical interactions of various telomeric structures on a short 
time scale and that is able to reproduce experimental measurements in mtDNA of yeast. Moreover, the model 
opens up a couple of interesting mathematical problems such as the validity of a quasi-steady state approximation 
and dynamic properties of discrete coagulation-fragmentation systems. We also identify and estimate key factors 
influencing the length distribution of telomeric circles, loops and strand invasions using numerical simulations.
}
{Received for publication X June 20XX and in final form X July 20XX.}
{Address reprint requests and inquiries to R.~Koll{\'a}r{, E-mail: kollar@fmph.uniba.sk.}}
\end{abstract}


\vspace*{2.7pt}
\begin{multicols}{2}
Telomeres are specialized nucleo-protein structures at the ends of linear DNA molecules involved in maintaining 
genomic stability \cite{mceachern2000telomeres}. Intensive research of telomeres started after Elisabeth Blackburn 
and Joseph Gall determined the sequences of termini of minichromosomes in the  macronucleus of the protozoan 
Tetrahymena \citep{blackburn1978tandemly}. Telomeres are associated with a number of proteins and RNAs playing 
various roles in essential functions of telomeres, such as: protection of DNA ends against degradation; masking 
the ends against inappropriate action of DNA repair machineries involved in processing of double-strand breaks; 
regulation of genes located in subterminal regions of chromosomes; localization of telomeric loci into specific 
regions of the nucleus; and pairing of homologous chromosomes during meiosis \cite{mceachern2000telomeres}.  

The observation that telomeric sequences of nuclear chromosomes in most eukaryotes consist of short \emph{tandem 
repeat units} (or shortly \emph{tandem repeats}) 
terminating with a $3'$ single-stranded overhang of the G-rich strand represented a basis for the 
discovery of mechanisms responsible for the solutions of the end-replication problem associated with the replication 
of linear DNA molecules catalyzed by conventional DNA polymerases \citep{olovnikov1971principle,watson2003origin}. 
The main mechanism responsible for replication of ends of nuclear chromosomes is based on \emph{telomerase}, the 
RNA-protein complex composed of the template RNA subunit  
%(\underline{te}lomerase \underline{R}NA; TER) 
and the protein catalytic subunit that extends the $3'$ single-stranded telomeric overhang and thus prevents 
shortening of chromosomes \citep{greider1985identification,greider1987telomere}. 


However, telomerase-mediated synthesis of chromosome ends is not the only mechanism of telomere 
maintenance \citep{lundblad2002telomere}. Examples of telomerase-independent pathways
\doiline

\noindent 
include (i) retrotransposition in {\it Drosophila} \citep{pardue2003retrotransposons}, 
(ii) telomeric loops (t-loops) \citep{griffith1999mammalian}, 
(ii) chromosome circularization in mutant strains of both fission yeast \citep{nakamura1998two} 
and Streptomyces \citep{qin2002survival} unable to maintain their telomeres and in mitochondria 
of {\it Williopsis} \citep{fukuhara1993linear} and {\it Candida} \citep{rycovska2004linear}, 
(iv) heterochromatin amplification-mediated telomere maintenance (HAATI) mechanism \citep{jain2010haati}, 
and (iv) homologous recombination. The latter was elaborated mainly by studies on telomere 
maintenance in yeasts \citep{teng1999telomere,teng2000telomerase}, but recently it was found 
to operate in a wide variety of organisms  including some insects, plants \citep{biessmann1997telomere} 
and humans \citep{dunham2000telomere}. Homologous recombination does not only help to maintain 
telomeres, but is also involved in generation of genomic plasticity and instability in the absence of telomerase, 
resulting in amplification and rapid deletion of telomeric DNA and  to the formation of extrachromosomal 
telomeric fragments \citep{lundblad2002telomere,niida2000telomere}. The importance of homologous 
recombination as means in telomere dynamics is underlined by the fact that it is one of the hallmarks of 
telomerase-deficient cancer cell lines maintaining their telomeres via \emph{Alternative Lenghtening of Telomeres}  (ALT)
mechanism \citep{cesare2010alternative,tomaska2009telomeric}.

The primary purpose of this study is a design and an analysis of a mathematical model of alternative telomere length maintenance 
that streamlines processes involved in ALT and thus offers a systematic view of the telomere interactions  
via an accurate description of biophysics involved in dynamics of telomere size distributions, i.e.,  factors that are often ignored in modeling. 
The model serves for investigations of parameters of ALT that are not measured by current experiments, and detects 
signatures of eventual individual biophysical components of telomere length maintenance.

{\it ALT mechanism in yeast mitochondria.}
A unique model system suitable for investigation of the replication and dynamics of telomeres maintained 
via ALT mechanism, and thus a source of experimental data for testing the mathematical model, 
occurs in mitochondrial DNA (mtDNA) of several yeast species including {\it C.~parapsilosis}, 
{\it C.~orthopsilosis}, {\it C.~metapsilosis}, {\it C.~salmanticensis} and {\it P.~philodendri} \citep{kovac1984yeast,nosek1995linear,nosek2004complete,rycovska2004linear}. In {\it C.~parapsilosis}, 
mitochondria contain two types of topologically distinct DNA replicons, linear and circular. Linear replicon 
corresponds to the mitochondrial genome carrying a standard set of mitochondrial genes and is represented 
by population of double-stranded DNA molecules (30,923 bp) terminating on both sides with inverted repeats 
consisting of a subterminal repeat (554 bp) and an array of tandem repeats (nx 738 bp). Individual molecules 
in this population differ in size by the number of full-length tandem repeats, ranging from 0 to more than 12. 
The shortest molecules (30,923 bp) terminate only with an incomplete tandem unit, while the termini 
of longer molecules (30,923 + nx 738 bp) are extended into arrays. The very ends of the molecules have 
single-stranded overhangs of about 110 bp, which in a fraction of molecules invade into the double-stranded 
region (\emph{invasion of the end}) thus forming a duplex loop structure termed t-loop \citep{tomaska2002t}. 

\begin{figurehere}
%\begin{figure}[htp] %R1
	\centering
	\includegraphics[width=9cm]{pic/tel_structures_3}
	\caption{Double stranded telomere replicons  \cite{tomaska2004alternatives}: 
$T_i$ a linear telomeric array with $i$ tandem repeats, and a 3' strand overhang, 
$C_j$ a telomeric circle with $j$ tandem repeats, 
$L_{i,j}$ a capped linear array (t-loop) formed by an invasion of the 3' end with $i+j$ tandem repeats 
out of which $j$ are part of the capping loop, 
and $Y_{i,j,k}$ a copy-choice structure (triple junctions) resembling a letter Y. 
\emph{Hotspots}--the particular locations within each tandem repeat most susceptible for recombinations are indicated.}
	\label{fig:CTL}
%\end{figure}
\end{figurehere}

The linear DNA molecules 
are accompanied by series of double-stranded circular molecules dubbed \emph{telomeric circles (t-circles)}. Their sizes 
correspond to integral multimers of the telomeric tandem repeat   with a length (nx 738 bp). The circular replicons may 
originate from the \emph{t-loop structures} and/or the \emph{telomeric arrays} at the ends of linear genomic molecules 
by recombination transactions followed by excision of a circle \citep{tomaska2000}. On the other hand,
the t-circles were shown to replicate independently of the linear genomic molecules via \emph{rolling-circle mechanism} 
leading to amplification of the linear arrays of tandem repeats. The amplified arrays can be resolved into t-circles. 
More importantly, they represent a substrate for recombinational mode of the mitochondrial telomere 
maintenance \citep{nosek2005amplification}. Investigation of mutant cells lacking the t-circles revealed that 
they contain a circularized derivative of the genome. This supports the idea that the t-circles play a key role 
in the mechanism of telomere maintenance \citep{kosa2006,rycovska2004linear}. This mechanism does not 
require telomerase activity; rather it relies on a relatively complex interplay among t-circles, lariat molecules 
(rolling-circle replication intermediates), t-loops and linear telomeric arrays either free (generated by rolling-circle mechanism) 
or attached to the termini of linear genomic molecules. In addition, a number of telomeric replicons
(i.e., telomeric tandem arrays, t-loops, t-circles, single-stranded overhangs) occur also at eukaryotic nuclear 
chromosomes pointing to a general significance of the mitochondrial model \citep{tomaska2009telomeric}. 
%%% add text about hotspots?


{\it Mathematical--biophysical model of ALT.}
The experimental and theoretical findings described above serve us as a basis for a construction of a predictive model of ALT dynamics.
Our key motivation arises from an aim to explain distinct properties appearing in size distributions of t-circles 
and linear telomeric arrays (see Fig.~\ref{fig:fit}). Species with relatively long tandem repeats {\it C.~parapsilosis} (nx 738 bp) \cite{} 
and  {\it P.~philodendri} (nx 288 bp) \cite{} seem to have exponentially decreasing 
distributions, whereas the size distribution for species  {\it C.~salmanticensis} with 
relatively short tandem repeats  (nx 104 bp) \cite{} is not monotone. 
Successful  reproduction of the most significant features of the experimental data without any  use of artificial 
free parameters suggests that the model captures fundamental elements of 
the ALT mechanism that are of a biophysical instead of a chemical or a biological nature. 
We claim that an interplay of entropy and elasticity is causing disparity in data with elastic forces decreasing 
a chance for existence and creation of short t-circles (with length less than approximately 500 bp) 
and entropy restricting the probability of emergence of long t-circles.  

%% scales
{\it Physical scales.}
An extreme complexity of biophysical and biochemical processes on various temporal and spatial scales involved in ALT creates a 
principal obstacle for modeling.  Although a model aiming to capture their full description, and particularly 
detailed molecular structure of telomeres and connecting proteins, may shed some light on individual 
aspects of telomere length maintenance, our approach is rather different. To achieve our primary 
goal to understand the size distributions of telomeric replicons we specify 
particular spatial and temporal scales of interest. We focus on a spatial scale with a typical size of hundreds 
of nanometers on which we distinguish size and overall shape (topology) of individual telomeres but ignore 
their detailed chemical structure, and particularly their nucleotide sequence. The typical time scale is set to a second, 
a relatively  fast time scale that is short for a DNA degradation, a cell replication, and for a significant contribution 
of synthesis of new tandem repeats via the rolling-circle mechanism.  
Although the fast time scale conservation of the total number of tandem repeats breaks down on  
a longer time scale commensurate with laboratory experiments, we argue that the measurements at each instance display 
``an equilibrium distribution"  for the given total number of tandem repeats, i.e. a quasi-steady state of the system. 


%% hints of results
{\it Overview of results.}
We construct a hierarchy of mathematical models of the alternative 
mechanism of telomere length maintenance with an increasing level of reduction and simplification. 
The robustness and reliability of our approach is demonstrated by a good agreement with experimental 
data  even on the simplest level that only accounts for interactions of t-circles. The compliance 
is reached for various yeast species in a model where the sole free parameter is the total number 
of tandem repeats in the system. The results of our numerical investigations of  the simplest model 
of the dynamics of the t-circle size distribution agree also with rigorous mathematical results 
both pointing out an interesting phenomenon of excess mass (gelation) possibly present in 
the biological system under certain specific conditions.  Due to their complexity, more realistic 
models that provide further information on size distributions of both linear telomeres and t-loops are 
only studied numerically. Our ultimate goal, modeling of ALT on a long time scale that accounts also 
for the rolling-circle mechanism, cell mitosis, and a DNA degradation, will be subject of our further experimental 
and mathematical research.



\section{Methods and Models} \label{methods}
The mathematical model of ALT mechanism is constructed successively in the following steps.  
First, we identify the key structures and reactions involved in ALT and 
determine a set of (biophysically and biologically) feasible assumptions about the experimental environment.
These postulates limit reaction types and enable us to determine a constitution of their kinetic rates that are built of
biophysical components for which, in the case of disagreement in the literature, we use a set of alternative formulae.  
Next, we perform reduction steps transforming the model to simpler forms better accessible to computational and theoretical investigations.  
Finally, a comparison with experimental data is used to select the appropriate alternatives for the individual components. 

{\it Assumptions.}
The key assumptions underlying the model are:
(i) all telomeres are double-stranded; 
(ii) more complex telomere replicons than those considered are dynamically unstable; 
(iii) telomeres do not have the super-helical structure (that changes their diffusive and reactive properties) and any interaction with such super-coiled structures is neglected; 
(iv) telomeres are well-mixed and are freely diffusing in a homogeneous three-dimensional confinement 
at a constant temperature (an assumption enabling a description of the dynamics 
in a form of system of ordinary differential equations rather than physically more realistic partial differential equations);
(v) all processes run on a short time scale on which synthesis of new tandem repeats and degradation of telomeres can be neglected and the total number of tandem repeats is conserved;
(vi) telomere interactions are local and site-specific;
(vii) there is an abundance of binding proteins in the environment, i.e., availability of proteins does not limit local reaction rates that can be then viewed 
as overcoming of energy barriers; 
(viii) there is only one specific reactive group (hotspot) per telomere tandem repeat that is susceptible to homological recombination
In general, more hotspots within one repeat can be present and their actual number can be treated as a parameter of the model 
that can be quantitatively estimated from the dynamical properties of the system. 
Obviously, in a real environment all of these assumptions fall short. Nevertheless, 
we assume that their failure does not significantly alter the dynamics of ALT, or
at least it does not change the qualitative behavior of the system. 

{\it Model set-up.}
In the absence of a fluid flow a motion of telomeres in an aquatic solution is governed 
by diffusion and elastic forces within the DNA. 
Despite the possibility that a tangled geometry of mitochondria may strongly influence diffusivity and bending of telomeres, 
small dimensions of telomeres (particularly of t-circles) and the experimental procedure that breaks internal cellular 
structure before each measurement allow us to neglect such spatial inhomogeneity. Furthermore, we describe the population 
dynamics of telomeres by chemical kinetics model expressed 
in a form of a system of ordinary differential equations \cite{phillips2009} that, in general,
requires large enough populations. However, in our application the
population levels of (large) telomeres with a large number of tandem repeats
are typically small, violating the condition. Nevertheless, we believe that
the good agreement of predictions of our model with experimental data provides its justification, 
although an extreme caution is necessary in interpretation of the results for 
populations of large telomeres.


{\it Reaction rates.}
The kinetic rates of telomere interactions (a bimolecular association, an unimolecular self-interaction, and dissociation)
can be structurally decomposed into individual factors common in the biophysical literature: 
\end{multicols}

\begin{itemize}
	\item Bimolecular association of free particles
		\begin{equation*}
		\text{RATE}_a = \begin{pmatrix}  \text{combinatorial} \\ \text{factor} \end{pmatrix} \times
				\begin{pmatrix}  \text{orientation} \\ \text{factor} \end{pmatrix} \times
				\begin{pmatrix}  \text{reaction-rate} \\ \text{limited factor} \end{pmatrix} \times
				\begin{pmatrix}  \text{diffusion-} \\ \text{limited rate} \end{pmatrix} 
		\end{equation*}
	\item Unimolecular self-interaction of a free particle (a reaction occurs between two parts of the same reactant)
		\begin{equation*}
		\text{RATE}_s = \begin{pmatrix}  \text{combinatorial} \\ \text{factor} \end{pmatrix} \times
				\begin{pmatrix}  \text{J-factor} \end{pmatrix} \times 2\times
				\begin{pmatrix}  \text{reaction-rate} \\  \text{limited factor} \end{pmatrix} \times
				\begin{pmatrix}  \text{diffusion-} \\ \text{limited rate} \end{pmatrix} 
		\end{equation*}
	\item Dissociation process (a reactant is decomposed into products)
		\begin{equation*}
		\text{RATE}_d = \begin{pmatrix}  \text{reaction-rate} \\  \text{limited factor} \end{pmatrix}\times
				\begin{pmatrix}  \text{diffusion-} \\ \text{limited rate} \end{pmatrix}  /
				{\begin{pmatrix}  \text{reaction}\\ \text{volume} \end{pmatrix}}
		\end{equation*}
\end{itemize}

\begin{multicols}{2}
The individual constituents are: 
(a) a combinatorial factor that accounts for the number of different configurations of reactants in 
a given reaction yielding the same set of products; 
(b) a diffusion-limited rate that measures the rate at which two perfectly spherical 
molecules encounter each other if driven only by a molecular diffusion; 
(c) a J-factor that accounts for bending 
of a telomere strand in terms of a local molar concentration of a reactive part at a given site; 
(d) an orientation factor that compensates for anisotropy of reactive surfaces of molecules; 
(e) a reaction-rate limited factor that accommodates reaction-limited interactions for homologous recombination and invasion of the end; 
and (f) a reaction volume equal to the sum of individual volumes of the reactants. 
Although for simplicity we neglect a presence of enzymes at reaction sites in the model, 
it is possible to build up these factors into our model, e.g., by modifying the orientation factor or a finite reaction rate.
For the detailed formulation of the model and its rates we refer the reader to the on-line supplement. 

{\it Size-dependent diffusion.}
Dependence of diffusion on polymer size \cite{robertson2006} is often suppressed in 
modeling of telomere interactions \cite{olofsson2010modeling,Peskin,rodriguez2010quantitative}
although it may represent an important factor in length maintenance, particularly for 
extra-chromosomal telomeric structures.
The classical theory of Flory \cite{FloryPolymer} suggests to neglect such an effect: {\it``%
%%the rate of reaction of a  functional group does not depend on the size of the molecule to which it is attached... 
While the range of diffusion of a terminal group within an interval of time which
is small compared to the period required for displacement of the molecule as a whole will be limited by its attachment to the polymer molecule, the group may nevertheless diffuse readily over a considerable region through rearrangements in the configurations of 
nearby segments of the chain."} 
Nevertheless, we consider the case of the size-dependent diffusion and compare it to 
to the case of the constant diffusion pointing out interesting phenomena worth a discussion. 
%Even in simplified models size-dependent diffusion improves convergence rates to equilibria 
%in the system and destroys 
%the metastable character of the dynamics. 
Note that a variable diffusion factor necessitates 
an inclusion of an orientation factor (specific site binding) \cite{Jackson} due to a disproportion in reactive-to-total-surface 
ratio of linear and circular telomere replicons.  Furthermore, presence of binding proteins and/or a variable molecular 
structure of telomeres require a consideration of reaction-rate limited kinetic rates \cite{phillips2009}.

{\it J-factor -- measure of DNA looping.}
The shape of telomeres is determined by random fluctuations (entropy) and elastic forces, 
influenced by their detailed molecular structure and binding proteins.  In biophysics all such 
effects are expressed in a quantifiable measure \emph{J-factor} \cite{Jacobson, EUreview} characterizing a local concentration 
of an interacting group relative to another one on a single DNA strand, enabling us
to determine individual reaction rates for self-interaction of DNA (recombination and self-end-invasion). 
The theory distinguishes between constrained and unconstrained (one of the reacting parts is at the free end of the polymer) 
looping of linear polymers, and (double) constrained looping of circular polymers. Unfortunately, in the large literature on J-factor 
there is no overall agreement on its values, particularly for constrained looping of circular polymers. 
In this letter we numerically test possible choices for J-factor and select one that best matches the experimental data. 
Thus our approach to data matching differs significantly from the traditional 
parameter fitting: instead of a parameter we fit a parameter-free model  with parameters from biophysical literature. 
This way our study may serve also a justification of selected models of J-factor. 

{\it Local interactions.}
The repetitive structure of telomeres conditions their possible mutual or self-interactions. We assume that they
can be classified into three types:  homologous recombination \cite{}, invasion of the end \cite{}, and dissociation (Fig.~\ref{fig:scheme}).
Hotspots (particular highly reactive parts of tandem repeats) are suspected to play a key role in a homologous recombination \cite{}.
Two identical recombination hotspots on the same or on two different strands under favorable conditions
(proximity and spatial alignment of reacting groups, availability of connecting proteins) 
form a transient complex resulting either in an original configuration of strands or in 
an exchange of  nucleotide sequences between the strands.

\begin{figurehere}
	\centering
	\includegraphics[width=8.3cm]{pic/scheme2_4}
	\caption{Schematic diagrams of local interactions: recombination, end invasion and dissociation of a copy-choice (triple junction) telomere replicons.}
	\label{fig:scheme}
\end{figurehere}

While homologous recombination is a reaction of two identical reacting groups (hotspots), 
the invasion of the end is an insertion of the 3' overhang at the terminus of a DNA strand into a non-terminal part 
of the same or another strand. Invasion into the same strand produces a capped form of a telomere (t-loop) while 
an invasion into a different strand yields a copy-choice structure (triple junction) of DNA strands. Again, the repetitive structure 
of the telomeres implies that each 3' overhang is only able to invade a particular location within a tandem 
repeat specific to the 3' overhang. Even under favorable conditions each local interaction of telomeres 
requires according to the Kramer's theory \cite{Jackson} crossing of a certain energy barrier (Fig.~\ref{fig_energy_barriers}). 
that effectively reduces the overall kinetic rate of reactions and consequently  introduces a finite reaction rate factor. 
Diffusion-limited rate of dissociation is calculated in a standard way \cite{Jackson} (see the on-line supplement for  more details). 
 \end{multicols}


\begin{figure}[htp] %1
	\centering
	\includegraphics[width = 15truecm]{pic/reactions4}
	\caption{List of chemical equations in the TCLY-model of ALT. The shadowed ovals represent a creation of a particular complex.}
	\label{fig_reactions}
\end{figure}

\begin{multicols}{2}


{\it CTLY-model.} 
Figure \ref{fig_reactions} schemtically displays interactions of four distinct types of telomere replicons: 
linear arrays, t-circles, t-loops, and triple junctions (see Fig.~\ref{fig:CTL}),
in our most complex model of alternative telomere length maintenance, 
dubbed the {\it CTLY-model}. The linear arrays are categorized according to their 
telomere length to the classes $T_i$, $i \ge 0$. Each linear array in $T_i$ terminates a linear mtDNA sequence 
and has $i$ full tandem repeats, with a free 3' single stranded overhang.
%%% typical for mtDNA as opposed to 5' overhangs typical for nDNA.
The (extra-chromosomal) circularized telomeres (t-circles)  are divided into classes $C_j$, $j \ge 1$, according to the number of 
full tandem repeats they contain.
Linear telomere replicons (capped or uncapped) belonging to four ends of a single chromosome are, for simplicity, 
considered as four independent particles. 
Although, as the reactions can produce larger and larger telomere structures, the size of the system is hypothetically unbounded 
(the system of reactions and various reactants of an arbitrary size is infinite), 
the initial value of the total number of telomeric tandem  repeats in the system creates an effective bound for the maximal size of a replicon.

Eventual telomere interactions in CTLY-model are limited to a homologous recombination of linear arrays and t-circles in reactions (A)-(C), 
and an invasion of an end of a linear array into the same or another telomere  in (D)-(F) (see Fig.~\ref{fig_reactions}). 
Two intermediate complexes are subsequently formed in (A)-(C) reflecting the three-stage process: 
approach of the two reactive parts, local protein mediated recombination, and dissociation.
We expect the intermediate complexes to have limited reactivity, so in the model they do not self-interact or recombine, 
and they are not invaded by telomeric overhangs to form more complex telomere replicons 
(double circles, quadruple junctions, etc.). 
%%This is consistent with the quasi-steady state assumption and will be discussed later. 



{\it Reduction to CT-model.}
It is common in mathematical models in biology and chemistry to reduce the size of the reacation system
by using a \emph{quasi-steady state approximation} (QSSA) \cite{Atkins, Eigen}. The QSSA can be used when a replacement of
a certain number of ordinary differential equations in the system (out of all)
by their equilibria does not significantly alter the dynamics of the whole system on a specificaly considered time scale. 
The components in the ``equilibrium" are often called \emph{slaved variables} as their approximated dynamics is slaved to the 
rest of the system by an algebraic law expressing the equilibrium in the reduced equations.  
Validity of such an approximation can be justified using the mathematical technique 
of singular perturbation \cite{SegelSlemrod,SM2003,TE2004} under a condition of a separation of time scales,  
specifically under overall imbalance of association and dissociation events. 
While the use of the QSSA can be more-or-less justified for every individual reaction 
on Fig.~\ref{fig_reactions}, the validity of the quasi-steady state needs 
to be justified for the full system of reactions simultaneously 
and not just for the individual reactions separately. Unfortunately, the existing QSSA 
theory does not cover such a case and thus it poses and interesting and difficult open 
mathematical problem. Nevertheless, here we assume QSSA in the CTLY-model as our numerical simulations suggest 
its validity for systems of large finite  sizes. 

The QSSA reduces the dynamics of the CTLY-model to a simplified 
\emph{CT-model}: 
\begin{eqnarray}
	 C_m + C_n  & \autorightleftharpoons{\small $k^{CC}_{f}$}{\small $k^{CC}_{b}$} &  C_{m+n}\, , \label{eq_CC}\\
	 T_m + C_n  & \autorightleftharpoons{\small $k^{TC}_{f}$}{\small $k^{TC}_{b}$} &  T_{m+n}\, , \label{eq_CT}\\
	 T_m + T_n  & \autorightarrow{\small $k^{TT}_{f}$}{} &  T_k + T_{m+n-k}\, ,\label{eq_TT}
\end{eqnarray}
where reaction (A) was reduced to (1),  (B), (D) and (E) were reduced and added up to \eq{eq_CT},
and similarly (C) and (F) were reduced and added up to \eq{eq_TT}. The final forward and backward rates 
and their derivation can be found in the supplement.


{\it Reduction to C-model.}
A closer look at the kinetic rates in \eq{eq_CC}--\eq{eq_TT} reveals a presence of considerably different 
time scales for fixed values of $m$ and $n$. Due to the size of the orientation 
factor that suppresses reactions of linear arrays t-circle interactions are the most favorable and thus 
also the most frequent as it is easier for the t-circles to orientate correctly. Therefore, one may organize the dynamics of  \eq{eq_CC}-\eq{eq_TT} according to the time scales in its indivudal reactions: \eq{eq_CT} determines the overall ratio of linear arrays and t-circles while \eq{eq_CC} 
controls the t-circle size distribution. Then \eq{eq_TT} plays the least important role, i.e., 
modifies the linear array distribution on even slower scale.

Therefore the distribution of the sizes of circular telomeres is dominated by the dynamics \eq{eq_CC} (that we call the {\it C-model}) that 
should be able to capture its essential features. The main advantage of reduction to C-model is that it fits the mathematical framework of  
coagulation-fragmentation processes \cite{AizBak,BakBak,BallCarr,daCosta,Fournier2004}, 
although the biophysically interesting reaction rates do not directly fall into categories studied theoretically up to this date. 
 \end{multicols}


\begin{figure}[tp] %R1
	\centering
	\includegraphics[width=12cm,height = 4cm]{pic/fit_all_modifDist2}
	\caption{The fit of the C-model to the experimentally measured t-circle distributions \cite{tomaska2000}.
	The fitted initial mass differs from the histogram apparent experimental mass by an amount not detected by the experiment.}
	\label{fig:fit}
\end{figure}

\begin{multicols}{2}
{\it Numerical Methods.}
The reaction kinetics equations in C-model and CT-model were represented by ordinary differential equations 
discretized using the Euler scheme (in numerical software MATLAB) that preserves the total number of tandem repeats in the system. 
The number of tandem repeats in any single linear array and t-circle in the implementation was limited to a finite size $2N$, where 
$N$ is sufficiently large to have minor influence on results creating an artificial boundary in the system. 
However, to avoid any significant effect of such a closure, 
only results for telomeres with less that $N$ full repetitions of the tandem repeat are considered to be reliable. 
Also, we interpret non-integer values of populations of linear arrays $T_i$ and t-circles $C_j$ as statistical measures of an average number
of particles of the given size and type in the sample.  
%The results of the numerical simulations were consequently compared to the biological experiments summarized in \cite{tomaska2000}.
A good agreement observed in a comparison of simulations of the full CTLY-model with the reduced CT-model 
demonstrates validity of our formal quasi-steady-state approximation.


\section{Results}

The main result of this essay is a confirmation that the proposed biological alternative mechanism of telomere length maintenance is in agreement with the experimental data as we were able to confirm both numerically and rigorously for models with a various level of simplification. Furthermore, our analysis uncovers several distinct phenomena that may be tested in the future biological experiements, and that may advance understanding of telomere length maintenance.

{\it Numerical results.}
The data displayed on Fig~\ref{fig:fit} confirm that the C-model is capable to capture the overall shape of the 
distribution of the t-circles for three distinct yeast species. After a selection of 
a particular reasonable values for physical paramaters and proper alternatives for modeling of J-factor that are the same for all three 
species (see the Section Parameters and Alternatives for details) the only fitted parameter in the numerical simulations 
is the total initial number of the tandem repeats that was not quantified in the experiment as the procedure 
does not guarantee that all the telomere replicons in the sample were detected. 
The relatively good agreement of a very rough model without fitting 
any artificial parameters also serves as a confirmation that the model captures the principal 
components of alternative mechanism of telomere length maintenance.
As expected the equilibrium distribution strongly depends on the species tandem repeat size due to the J-factor influence. 
Shorter telomeres have a greater difficulty to bend (due to elastic forces) while longer (pliable) strands have smaller tendency 
for circularization due to a small probability of a proper alignment of the specific reactive sites (entropy effects). 
Due to an algebraic homogeneity of a corresponding factor, the equilibrium distributions for the C-model using 
the  size-dependent and constant diffusion are the same, and the only difference between them in dynamics appears only in 
the rate of convergence of an initial distribution to the equilibrium. 

Numerical simulations with large initial number of tandem repeats demonstrate a presence of the phenomenon 
of {\it excess mass} (gelation). In cases when the initial condition prescribes a number of tandem repeats 
that cannot be accomodated by an equilibrium of the system of  the highest possible mass,  the extra 
(excess) mass is transfered by the dynamics of the system to higher and higher modes, i.e.,  increasing 
a probability of presence larger and larger t-cirlces, while the distribution of telomeres of fine sizes converges 
to the equilibrium of the largest possible mass. Such a case appears to happen any time the initial  total count 
of tandem repeats iexceeds a critical value corresponding to a dynamical equilibrium of the highest possible mass.
This mechanism may be responsible for creation of the extra long t-tcircles in situations when the internal clock
 of the cell breaks down, the cell produces more and more telomere tandem repeats uncontrollably by the rolling circle mechanism, and 
an extremely large t-circle is created dynamically in the sample. Whether this phenomenon correspond to the unusually large experimentally observed t-circles (e.g. a t-circle with 9 full tandem repeats in {\it C. parapsilosis} \cite{}) is unclear and remains a subject of our further investigations.

Furthermore, our numerical simulations of the C-model show a significant speed up of the rate of convergence 
to the equilibrium for systems with size-dependent diffusion compared 
to the systems with constant diffusion (Fig.~\ref{fig:conv}). That may be a sign of
the fact that size-dependent diffusion in the biological application enhances adaption 
to changes of environment and other types of stimuli. Such a difference in dynamics can serve as a lead
to an answer to the question whether  the size-dependent diffusion corresponds has a physical significance for 
the behavior of the system or it can be neglected as the theory of Flory \cite{FloryPolymer} suggests. 

{\it Analytical results.}
Convergence to equilibrium in C-model for kinetic rates for both size-dependent and constant
diffusion was also established by rigorous mathematical analysis \cite{Carr}. 
For small initial mass
both types of kinetic rates lead to the same exponentially decaying (up to an algebraic factor) population 
distribution of t-circles. In both cases the equilibrium distribution is a multiple of the J-factor that justifies our hypothesis that the J-factor is the primary reason for the different overall shape of the experimentally measured size distributions of t-circles for various species of yeast.  

The rigorous analysis also reveals the critical value for the initial number of tandem repeats that corresponds 
to an equilibrium of the highest possible mass---the unique equilibrium with an algebraic decay, 
and thus a threshold for the phenomenon of excess mass.  The analysis confirms numericaly observed behavior: 
the distribution of telomeres of fine sizes (weakly) converges to the equilibrium of the largest possible mass 
and the extra mass moves to larger and larger t-circles, eventually leaving the system through a t-circle of ``an infinite size".  
While such a rigorous result can be established for the system with constant diffusion, the size-dependent diffusion 
creates a significant obstacle to existing methods of analysis and poses a challenging theoretical mathematical problem, 
although the numerical simulations suggest that the results for both diffusion alternatives should be the same. 
\end{multicols}

\begin{figurehere} %R1
	\centering
	\includegraphics[width=18cm]{pic/beta_gamma}
	\caption{(left panel) The size distribution dependence on parameter $\beta$. 
	The blue line corresponds to value of $\beta\rightarrow \infty$ whereas the orange line corresponds to $\beta\rightarrow 0$. 
	(right panels) The size distribution dependence on the size of parameter $\gamma$ with $\gamma\rightarrow 0$ (green line) 
	and $\gamma\rightarrow \infty$ (red line). Initial conditions are chosen so that the initial mass in t-circles and linear 
	telomeres are both equal to the equilibrium mass.}
	\label{fig:res5}
\end{figurehere}


\begin{figurehere}
\centering
\includegraphics[width=11cm]{pic/convergence}
\caption{A comparison of a convegernce of a t-circle distributionto an equilibrium for systems with size-dependent (Case 2) and constant diffusion (Case 3) on a linear (left panel) and a logarithmic (right panel) scale. In both cases the same initial data with the total number of telomeric repeats less than critical is used. The critical solution corresponds to an equilibrium with a maximum possible mass while the exact solution is the exact equilibrium of the system with th total number of tandem repeats corresponding to the initial data, in both cases the exact solutions is reached as $t \rightarrow \infty$.}
\label{fig:conv}
\end{figurehere}


\begin{multicols}{2}
\section{Parameters and Alternatives}

{\it Non-dimensional parameters.}
In addition to known biophysical parameters that can be found in the literature the kinetic rates in the
C-model and TC-model depend on three non-dimensional local molecular parameters $\alpha, \beta$ and $\gamma$. 
They depend on a relative size of the correction factors to the diffusion-limited kinetic rates $D$ in the individual reactions 
depicted on Fig.~\ref{fig_energy_barriers} that correspond to energy barriers needed to form or break a bond in 
a recombination ($g$) and an end invasion ($g^\ast$) with the appropriate lower indices indicating the type of reaction. 


\begin{figurehere}
	\centering
	\includegraphics[width = 8.3truecm]{pic/energy_barriers_2}
	\caption{Kinetic rates of local reactions in the end invasion and recombination.}
	\label{fig_energy_barriers}
\end{figurehere}

The parameter $\alpha = g_1^{\ast}/g_1$ is the ratio of 
correction terms to diffusion-limited kinetic rates of invasion of the end and homologous recombination measured 
in terms of the difference in height of the corresponding energy barriers. If $\alpha < 1$ recombination is 
an energetically preferred mechanism to the end invasion, while $\alpha > 1$ otherwise.  
The lack of experimental measurements of $\alpha$ motivates our numerical parametric study confirming 
that $\alpha$ has a very minor influence on  the dynamics of the system.  
We understand that this is a consequence of the fact that each telomere 
has only one end but possibly many hotspots (one in each tandem repeat) and thus 
an invasion is, in general,  a rare event compared to recombination. We set $\alpha =e^{0.3}$ based on 
\cite{}.  

The ratio of correction factors to diffusion-limited kinetic rates that measures 
the difference in height of energy barriers in recombination 
and dissociation of recombination complex is denoted as $\beta = g_{-1}/g_0$. 
The value of $\beta$ reflects the relative speed of recombination complex dissociation to the speed 
of recombination once the complex is formed. We assume a very fast recombination ($\beta=0$) in the remaining 
numerical simulations but here we study the effects of $\beta$ for all admissible values $\beta\in(0,\infty)$. 

The numerical results on Figure~\ref{fig:res5} show that the size distribution of t-circles may be significantly affected 
by $\beta$ for {\it C.~salmanticensis}.  unlike the linear arrays distribution that remains almost unchanged (not shown). 
The yeast with a longer tandem repeat ({\it C.~parapsilosis} and {\it P.~philodendri}) show almost no sensitivity with respect 
to $\beta$ that may be explained by a better stability of the recombination complex due to less significant elasticity effects. 
The minor changes in the size distributions involve a slight shift of the distribution accompanied by the slight change of 
the exponential decay. Since we do not have accurate experimental data to measure the rate of this exponential decay 
we do not study these minor differences.

Finally, $\gamma = g^{\ast}_{-1}/g^{\ast}_2$ is the ratio of the correction factors to diffusion-limited rates measuring the difference in 
energy barriers of breaking of a strand invasion in a two possible ways. Values $\gamma \ll 1$ 
represent a forward reaction where the invasion of the overhang is followed by a break of the invaded strand. 
In the case of a linear telomere invading to a t-circle this results into telomere lengthening. However, in the case
of a self-invasion this results into telomere shortening. Therefore it is not clear if the change of gamma results in 
a mass redistribution into t-circles or into linear telomeres. The numerical results shown on Figure~\ref{fig:res5} 
show that a decrease of $\gamma$ results into more mass in linear telomeres and less in t-circles. At the same time 
the mean of the t-circle size distribution shifts to lower values in accordance with the results in part (a) for a simpler C-model.

On the other hand, $\gamma \gg 1$ represents a dominating reversible reaction for the end invasion. Numerical simulations indicate that with an increasing $\gamma$ the total amount of tandem repeats in t-circles increases. As a consequence, the mean of the size distribution will shift  to larger values. We observed no significant change in the t-circle or linear telomere size distributions for other studied 
yeast organisms ({\it C.~parapsilosis} and {\it P.~philodendri}) with respect to the parameter $\gamma$. 
We set $\gamma = 1$.


\noindent
{\it Dimensional biophysical parameters.}
There are two parameters that are specific for individual speciec of yeast and that have known values: 
tandem repeat length that was discussed in details above and the total chromosome length that influences diffusive properties of linear arrays in size-dependent diffusion alternative. The magnitude of the local correction do diffusion-limited kinetic rate $g_1$ (see Fig.~\ref{fig_energy_barriers}) only affects the overall time scale of the system  and does not influence the system dynamics, so by a proper rescaling of time we eliminate it from the system. 

{\it Modeling alternatives.}
The numerical experiments were performed for three yeast species: 
{\it C.~parapsilosis}, {\it P.~philodendri} and {\it C.~salmanticensis}, differing in the total genome size 
and, more imporantly,  in the length of the tandem repeat.  
The J-factor is modeled using three alternatives described in the supplement, 
namely: (I) separate models for constrained and unconstrained J-factor \cite{Shimada}; 
(II) a universal model {\it in vitro} \cite{Han}; and (III) and a universal model {\it in vivo},  \cite{Han}. 
Additionally, we consider three alternatives for the circular J-factor: (A) Gaussian model with detailed balance 
using a linear J-factor; (B) the model of \cite{Rippe1998} 
with $d=0$ and (C) with $d=0.13$, where $d$ measures the length of a mediating protein in nm. 
Note that this creates in total 9 options for combination of linear and circular J-factor in the TC-model, although 
it reduces to 5 options in C-model, as the linear $J$-factor occurs there only in the Gaussian alternative. 

There were two alternatives for the diffusion coefficients in the kinetic rates, one with constant diffusion 
(corresponding to the expected diffusion coefficient of a smallest possible t-circle) and another 
one with size-dependent diffusion. For completeness, we have also made a comparision of the dynamics with the so-called 
Becker-D{\"o}ring model where an interaction of particles requires that one of the reacting species is of the size one (tandem repeat). 

To capture the overall shape of the t-circle distribution we compared experimental data with the computer simulations 
for the C-model using different J-factor alternatives. The comparison was done for yeast with short tandem repeats 
({\it P.~philodendri} and {\it C.~salmanticensis}) where the J-factor forms differ (see Figure~\ref{fig:res1}, right panel), 
while keeping the total tandem repeat count mass the same in all simulations. 

The  results for {\it C.salmanticensis} on Figure~\ref{fig:res1} yield a good match for two models: the Gaussian circular J-factor 
combined with an {\it in vivo} linear J-factor, model III; and a Rippe formula with $d=0.13$. Other equilibria incorrectly 
predict very small density of the shortest t-circles.

The results for both yeast species imply that the Rippe formula that takes into account existence of a mediating protein matches 
the t-circle distribution better. Furthermore, we performed additional numerical study on the CT-model distinguishing 
nine J-factor alternatives (results on Figure~\ref{fig:res15}) confirming the optimality of {\it in vivo} linear J-factor 
for all circular J-factors and the best fit being the combination of a Rippe formula ($d=0.13$) with the {\it in vivo} 
J-factor (model III). 

\end{multicols}

\begin{figurehere}
	\centering
	\includegraphics[width=17cm]{pic/Cplot25_fitJ_modifDist}
	\caption{A comparison of numerical simulations of the C-model with the experimental data \cite{tomaska2000} for various modeling alternatives. The total number of the tandem repeats is the same as detected in the experiment.}
	\label{fig:res1}
\end{figurehere}


\begin{figurehere}
	\centering
	\includegraphics[width=13cm]{pic/CShist33_fitJ_modifDist}
	\caption{A comparison of numerical simulations of the CT-model with the experimental data \cite{tomaska2000} for various modeling alternatives. Different choices for the, respectively, linear and circular J-factor are in columns and rows, respectively. 
The total number of the tandem repeats is the same as detected by the experiment.}
	\label{fig:res15}
\end{figurehere}

\vspace{\baselineskip}

\begin{figurehere}
	\centering
	\includegraphics[width=16cm]{pic/C_InitialMass}
	\caption{A dependence of the equilibrium t-circle distribution in the C-model on the total number of tandem repeats
	 (with the logaritmic plot on the right panel).  Data are renormalized to a total initial mass equal to one.}
	\label{fig:res2}
\end{figurehere}

\vspace{\baselineskip}

\begin{figurehere}
\centering
\includegraphics[width=0.4\textwidth]{pic/ratio}
\includegraphics[width=0.4\textwidth]{pic/Tnumber}
\caption{A dependence of the fraction of t-circles (left panel) and the shape of the t-circle distribution (right panel) 
in the TC-model on the total number of linear arrays. With an increasing number of linear arrays the ratio of t-circles decreases exponentially.  The shape of the t-circle distribution in the TC-model does not vary significantly with an increasing total number of linear arrays. }
\label{fig:res3}
\end{figurehere}


\begin{multicols}{2}
{\it Dependence of equilibria on experiment specific parameters.}
There are two additional experiment specific parameters that are conserved under the dynamics: total number of tandem repeats and the total number of linear arrays in the experiment. They both significantly influence the equilibrium distributions of t-circles (in C-model) and both t-circles and linear arrays (in CT-model). 

The parameter of initial mass corresponds to the size of the experiment. However, only a fraction 
of the t-circles in the experiment were detected and thus the rate of detection is an unknown constant of the experiment, 
related to the total mass and the detected mass via a linear relationship.  

Because the C-model is nonlinear with respect to the t-circle concentrations the stationary density should depend on the 
size of the initial mass in the system. This can be observed on Figure~\ref{fig:res2} where the stationary solutions rescaled 
to a unit total mass are compared. The larger initial mass implies a shifted mean of the distribution towards larger structures. 
Moreover, the tail of the distribution is exponential which is consistent with the theoretical results of the 
coagulation-fragmentation problems in slightly different regimes (studied elsewhere).

Tandem repeats rearrange into circular and linear replicons in such a way that both the distribution of t-circles and linear arrays
show exponential tails. Moreover, the large number of linear telomeres allows  more material to be distributed into 
the ends of linear telomeres as opposed to t-circles thus changing the relative mass in t-circles, see Figure~\ref{fig:res3}. 
However, not only the relative telomeric mass changes with the telomere count but also the shape of the distribution 
changes slightly due to the nonlinear nature of the problem. The ratio $C/(C+T)$ in equilibrium can be well approximated 
by an exponential $\exp(-cT)$ for some constant $c>0$. 

%\begin{figurehere}
%\centering
%\includegraphics[width=0.4\textwidth]{pic/ratio}
%\caption{Dependence of the fraction of t-circles in the TC-model on the total number of linear arrays. With an increasing number of linear 
%arrays the ratio of t-circles decreases exponentially. }
%\label{fig:res3a}
%\end{figurehere}
%
%\begin{figurehere}
%
%	\centering
%	\includegraphics[width=0.4\textwidth]{pic/Tnumber}
%	\caption{The shape of the  t-circle distribution in the TC-model does not vary significantly with an increasing total number of linear arrays. }
%	\label{fig:res3}
%\end{figurehere}


\section{Discussion}
The rigorous analysis and numerical simulations reveal an interesting feature possibly exploited 
in estimation of quantifiable detection rate of the PCR method. Indeed, the equilibrium of the system 
has a nonlinear response to the number of tandem repeats in the sample. If one interprets the experimentally 
obtained southern blot samples as relative quantities of measured species of different sizes, the model then 
identifies the total number of  tandem repeats in the sample by using the location and a relative size of the peak 
of the distributions of linear arrays and t-circles. It allows us to relate the results of the PCR method to an approximate number of particular polymers in the sample. Clearly, such a method would not be possible if the system response would be linear. 

Extensive studies in biophysical literature identify a dependence of the individual factors of the kinetic rates 
on various environmental parameters: temperature, salinity, or pH. Variations of these environmental factors in experiments and the subsequent measurement of the telomere size distributions, together with the 
comparison with the predictions of the model, can eventually serve as a further justification of the 
individual module, and also can lead to the further factor-dependent enhancement of the model. 

One of the key issues of our study in an assumption on homogeneous environment in the experimental sample. In fact, the 
measurement is taken from a collection of cells and thus the resulting distribution can be viewed as a sum of individual distributions 
in the cells. On the other hand, all the telomeres in the sample coexist together for a short time after a break of 
cellular structures ({\it in vitro}) in a solution where they can freely diffuse.  To better understand the telomere length 
maintenance process in a single cell rather than in a collection of cells (with distributed phases of cell cycle), single cell measurements must be obtained experimentally. 

Due to a lack of available dynamical experimental data (time series) the major assumption of this study is an assumption on a time scale. 
We have only considered a very fast time scale on which the total number of tandem repeats in constant. 
Nevertheless, our pricipal goal is to mode the alternative mechanism of telomere maintenance as a whole on a longer 
time scale that includes sythesis of tandem repeats via rolling-circle mechanism, cell mitosis, and subsequently 
the end-replication problem, and the DNA degradation. Such a complex model may help to answer many imporatant biological 
question as how many of tandem repeats are transfered from a mother cell to a dauthger cell and in what form, how the internal 
clock of a cell influence the telomere maintenance, etc., but it necessitates an availability data from new experiments. 


\section*{SUPPLEMENTARY MATERIAL}

An on-line supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.\vspace*{6pt}

\bibliography{ref_telomeres}

\end{multicols}

\end{document}

